Read in the packages. The working directory is wherever the R Notebook is located.
Read in the spreadsheet and take a look at the data.
###read in spreadsheet
loc <- read_xlsx("Original-spreadsheets/all species New_6-14-19.xlsx") %>%
janitor::clean_names() %>%
mutate(reproductive_mode = as.factor(reproductive_mode))
#get the number of individuals, and the sexuality counts per species
count_repro_mode <- loc %>%
group_by(genus, species, reproductive_mode) %>%
dplyr::count() %>%
mutate(genus_species = str_c(genus, species, sep = "_"),
genus_species = str_replace_all(genus_species, " ", "_"),
genus_species = str_replace_all(genus_species, "\\.", "")) %>%
ungroup() %>%
mutate(genus_species = fct_reorder(genus_species, n, sum)) %>%
ggplot(aes(x = genus_species, y = n, fill = reproductive_mode)) +
geom_col() +
coord_flip() +
theme_minimal()
count_repro_mode##Map Plot a leaflet map of the localities. The leaflet map is interactive. You can click on the localities and a flag with some metadata will pop up!
#make locality shape file and assign WGS coord system
coord_points <- st_as_sf(loc, coords = c("longitude", "latitude"),
crs = 4326, agr = "constant")
#use sourced plot_locs_leaflet script to plot localities
all_plot <- plot_locs_leaflet(loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(all_plot, url = paste0(getwd(), "/plots/repro_mode_plots/all_species_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/all_species_map.pdf"))Obtain the bioclim layers for analysis. I’m using all 19 for this preliminary exploration. I plotted the first bioclim just to make sure nothing seems wonky. I’m using CHELSA data downloaded from their website. Since the files are huge, I only unzip them one at a time, crop them, and write them to GeoTiff files that I can then load in as a rasterstack.
##get chelsa data
#chelsa_folder <- "/Users/connorfrench/Dropbox/Old_Mac/climate-data/chelsa_30s_bio"
#zip_files <- list.files(chelsa_folder, full.names = TRUE)
#using the Unarchiver commandline tools for Mac to unzip the 7zip chelsa layers. Regular unzip() does not work with 7z zipped files
#for (file in zip_files) {
#set temp directory
# tempd <- tempdir()
# system(paste("unar", file, "-o", tempd))
# r <- raster(list.files(tempd, pattern = "*.tif", full.names = TRUE)) %>%
# crop(extent(166, 179, -48, -34))
# writeRaster(r, filename = paste0("~/Desktop/", list.files(tempd, pattern = "*.tif")), format = "GTiff")
# unlink(tempd, recursive = TRUE)
#}
clim_files <- "/Users/connorfrench/Dropbox/Old_Mac/climate-data/chelsa_30sec_NewZealand/chelsa_bioclims_NZ.tif"
w <- stack(clim_files)This is a PCA of the climate data extracted for each locality, rather than a PCA of the total climate space.
Run the pca and check out variable loadings and proportion of variance explained by components.
#extract data from worldclim for each locality. Making this into a data frame with columns labeled so the row labeling lines up after I remove the NAs.
#extract data from worldclim for each locality.
coords <- data.frame(latitude = loc$longitude, longitude = loc$latitude)
loc.clim <- dplyr::bind_cols(loc, raster::extract(w, coords, method = "simple", df = TRUE)) %>%
drop_na(chelsa_bioclims_NZ.1) %>%
dplyr::select(-ID)
#make a matrix of only bioclim values
clim.mat <- loc.clim[,grep("bio", names(loc.clim))] %>% as.matrix()
#run pca on climate variables
clim.pca <- prcomp(clim.mat, scale = TRUE)
summary.pca <- summary(clim.pca) #check out the components
#plot tables
summary.pcaImportance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7
Standard deviation 2.9511 2.4947 1.36934 1.02403 0.80149 0.54359 0.32182
Proportion of Variance 0.4584 0.3276 0.09869 0.05519 0.03381 0.01555 0.00545
Cumulative Proportion 0.4584 0.7859 0.88462 0.93981 0.97362 0.98917 0.99462
PC8 PC9 PC10 PC11 PC12 PC13
Standard deviation 0.24972 0.11677 0.10148 0.08340 0.05671 0.04746
Proportion of Variance 0.00328 0.00072 0.00054 0.00037 0.00017 0.00012
Cumulative Proportion 0.99790 0.99862 0.99916 0.99953 0.99970 0.99982
PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.03910 0.03147 0.02392 0.01408 0.01153 0.00911
Proportion of Variance 0.00008 0.00005 0.00003 0.00001 0.00001 0.00000
Cumulative Proportion 0.99990 0.99995 0.99998 0.99999 1.00000 1.00000
| PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | -0.301 | 0.148 | 0.174 |
| chelsa_bioclims_NZ.2 | -0.285 | 0.143 | 0.276 |
| chelsa_bioclims_NZ.3 | -0.310 | 0.142 | 0.097 |
| chelsa_bioclims_NZ.4 | -0.150 | -0.358 | 0.010 |
| chelsa_bioclims_NZ.5 | -0.189 | -0.327 | 0.024 |
| chelsa_bioclims_NZ.6 | -0.135 | -0.362 | -0.003 |
| chelsa_bioclims_NZ.7 | -0.169 | 0.136 | 0.089 |
| chelsa_bioclims_NZ.8 | -0.188 | -0.328 | 0.031 |
| chelsa_bioclims_NZ.9 | -0.130 | -0.365 | -0.006 |
| chelsa_bioclims_NZ.10 | -0.119 | -0.361 | -0.036 |
| chelsa_bioclims_NZ.11 | -0.214 | -0.280 | 0.094 |
| chelsa_bioclims_NZ.12 | 0.277 | -0.106 | 0.360 |
| chelsa_bioclims_NZ.13 | 0.261 | -0.105 | 0.306 |
| chelsa_bioclims_NZ.14 | 0.248 | -0.083 | 0.383 |
| chelsa_bioclims_NZ.15 | -0.219 | 0.123 | 0.492 |
| chelsa_bioclims_NZ.16 | -0.314 | 0.143 | 0.009 |
| chelsa_bioclims_NZ.17 | 0.272 | -0.106 | 0.380 |
| chelsa_bioclims_NZ.18 | -0.167 | 0.128 | -0.084 |
| chelsa_bioclims_NZ.19 | -0.247 | 0.100 | 0.317 |
Two plots: One plot of the PCA colored according to genus, with convex hulls surrounding the genera. It looks like this reflects a latitudinal gradient in temperature! You can interact with the PCA plot by clicking on points to view associated metadata. You can isolate the genus you want to view by double clicking the genus in the legend! You can also remove a genus by clicking on it once. There’s some other functionality you can explore in the toolbar at the top of the plot. The second plot is a PCA colored according to reproductive mode. It looks like asexual populations occupy slightly larger niche space, but both reproductive modes have a similar niche center.
#add pca results to loc.clim data frame
loc.clim <- data.frame(loc.clim, clim.pca$x)
#use sourced plot_clim_pca function to plot the pca results. args are the data set with species names and PC axis values and the pca summary
all_pca <- plot_clim_pca(loc.clim, summary.pca, factor = "genus")Ignoring unknown aesthetics: text
#use sourced plot_clim_pca function to plot the pca results. args are the data set with species names and PC axis values and the pca summary
repro_pca <- plot_clim_pca(loc.clim, summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#save the plot colored by genus
#htmlwidgets::saveWidget(all_pca, paste0(getwd(), "/plots/repro_mode_plots/all_species_pca_genus.html"), selfcontained = TRUE)
#save the plot colored by reproductive mode
#htmlwidgets::saveWidget(repro_pca, paste0(getwd(), "/plots/repro_mode_plots/all_species_pca_repro.html"), selfcontained = TRUE)These are PCAs of environmental space for species within genera. Each climate PCA is of localities for a single genus, colored by species. I’m doing this even for genera with one species, so it’s easy to see if certain localities group together.
#source function to conduct a PCA on individual species
summary.list.acan <- species_pca_fun(loc.clim, "acanthoxyla")
#plot
acan_plot <- plot_clim_pca(summary.list.acan$loc.clim, summary.list.acan$summary.pca, "reproductive_mode")Ignoring unknown aesthetics: text
#save pca plot
#htmlwidgets::saveWidget(acan_plot, paste0(getwd(), "/plots/repro_mode_plots/acanthoxyla_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
acan_loc <- loc %>%
filter(genus == "acanthoxyla")
#use sourced plot_locs_leaflet script to plot localities
acan_map <- plot_locs_leaflet(acan_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(acan_map, url = paste0(getwd(), "/plots/repro_mode_plots/acan_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/acan_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7
Standard deviation 2.8059 2.5135 1.7771 0.88976 0.66876 0.44564 0.36423
Proportion of Variance 0.4144 0.3325 0.1662 0.04167 0.02354 0.01045 0.00698
Cumulative Proportion 0.4144 0.7469 0.9131 0.95475 0.97829 0.98874 0.99573
PC8 PC9 PC10 PC11 PC12 PC13
Standard deviation 0.22139 0.1064 0.09506 0.06761 0.05169 0.03848
Proportion of Variance 0.00258 0.0006 0.00048 0.00024 0.00014 0.00008
Cumulative Proportion 0.99831 0.9989 0.99938 0.99962 0.99976 0.99984
PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.03336 0.02780 0.02371 0.01891 0.01400 0.01128
Proportion of Variance 0.00006 0.00004 0.00003 0.00002 0.00001 0.00001
Cumulative Proportion 0.99989 0.99993 0.99996 0.99998 0.99999 1.00000
loadings.acan <- summary.list.acan$summary.pca$rotation
knitr::kable(round(loadings.acan[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | -0.266 | -0.230 | 0.170 |
| chelsa_bioclims_NZ.2 | -0.238 | -0.200 | 0.297 |
| chelsa_bioclims_NZ.3 | -0.277 | -0.240 | 0.056 |
| chelsa_bioclims_NZ.4 | 0.262 | -0.267 | 0.039 |
| chelsa_bioclims_NZ.5 | 0.236 | -0.288 | 0.068 |
| chelsa_bioclims_NZ.6 | 0.268 | -0.258 | 0.015 |
| chelsa_bioclims_NZ.7 | -0.192 | -0.051 | 0.217 |
| chelsa_bioclims_NZ.8 | 0.239 | -0.286 | 0.073 |
| chelsa_bioclims_NZ.9 | 0.269 | -0.256 | 0.012 |
| chelsa_bioclims_NZ.10 | 0.272 | -0.249 | -0.011 |
| chelsa_bioclims_NZ.11 | 0.191 | -0.283 | 0.148 |
| chelsa_bioclims_NZ.12 | 0.219 | 0.232 | 0.285 |
| chelsa_bioclims_NZ.13 | 0.232 | 0.246 | 0.151 |
| chelsa_bioclims_NZ.14 | 0.128 | 0.138 | 0.463 |
| chelsa_bioclims_NZ.15 | -0.150 | -0.097 | 0.477 |
| chelsa_bioclims_NZ.16 | -0.278 | -0.245 | -0.034 |
| chelsa_bioclims_NZ.17 | 0.197 | 0.197 | 0.370 |
| chelsa_bioclims_NZ.18 | -0.074 | -0.304 | -0.109 |
| chelsa_bioclims_NZ.19 | -0.246 | -0.049 | 0.321 |
#conduct pca
summary.list.argo <- species_pca_fun(loc.clim, "argosarchus")
#plot
argo_plot <- plot_clim_pca(summary.list.argo$loc.clim, summary.list.argo$summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. It's an annoying bug that hasn't been fixed yet.
#htmlwidgets::saveWidget(argo_plot, paste0(getwd(), "/plots/repro_mode_plots/argosarchus_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
argo_loc <- loc %>%
filter(genus == "argosarchus")
#use sourced plot_locs_leaflet script to plot localities
argo_map <- plot_locs_leaflet(argo_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(argo_map, url = paste0(getwd(), "/plots/repro_mode_plots/argo_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/argo_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7
Standard deviation 2.8508 2.4341 1.6389 1.09192 0.76813 0.51507 0.35817
Proportion of Variance 0.4277 0.3118 0.1414 0.06275 0.03105 0.01396 0.00675
Cumulative Proportion 0.4277 0.7396 0.8809 0.94367 0.97473 0.98869 0.99544
PC8 PC9 PC10 PC11 PC12 PC13
Standard deviation 0.2222 0.10883 0.10035 0.08126 0.05441 0.04318
Proportion of Variance 0.0026 0.00062 0.00053 0.00035 0.00016 0.00010
Cumulative Proportion 0.9980 0.99866 0.99919 0.99954 0.99970 0.99980
PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.03895 0.03285 0.02417 0.01722 0.01643 0.01134
Proportion of Variance 0.00008 0.00006 0.00003 0.00002 0.00001 0.00001
Cumulative Proportion 0.99988 0.99993 0.99996 0.99998 0.99999 1.00000
loadings.argo <- summary.list.argo$summary.pca$rotation
knitr::kable(round(loadings.argo[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | -0.266 | -0.227 | -0.182 |
| chelsa_bioclims_NZ.2 | -0.241 | -0.208 | -0.293 |
| chelsa_bioclims_NZ.3 | -0.279 | -0.228 | -0.110 |
| chelsa_bioclims_NZ.4 | -0.225 | 0.312 | -0.037 |
| chelsa_bioclims_NZ.5 | -0.257 | 0.272 | -0.023 |
| chelsa_bioclims_NZ.6 | -0.202 | 0.330 | -0.026 |
| chelsa_bioclims_NZ.7 | -0.141 | -0.175 | 0.021 |
| chelsa_bioclims_NZ.8 | -0.255 | 0.275 | -0.033 |
| chelsa_bioclims_NZ.9 | -0.197 | 0.334 | -0.030 |
| chelsa_bioclims_NZ.10 | -0.176 | 0.347 | 0.013 |
| chelsa_bioclims_NZ.11 | -0.277 | 0.219 | -0.090 |
| chelsa_bioclims_NZ.12 | 0.266 | 0.149 | -0.318 |
| chelsa_bioclims_NZ.13 | 0.267 | 0.108 | -0.280 |
| chelsa_bioclims_NZ.14 | 0.227 | 0.162 | -0.310 |
| chelsa_bioclims_NZ.15 | -0.113 | -0.143 | -0.516 |
| chelsa_bioclims_NZ.16 | -0.291 | -0.224 | -0.009 |
| chelsa_bioclims_NZ.17 | 0.258 | 0.163 | -0.324 |
| chelsa_bioclims_NZ.18 | -0.098 | -0.073 | 0.247 |
| chelsa_bioclims_NZ.19 | -0.191 | -0.155 | -0.384 |
Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.
#get background env't for the species
ahor_bg_env <- bg_env_crop(argo_loc,
species = "horridus",
environment = w,
buffer = 0.5)
#enfa for the sexual species
ahor_sexual_enfa <- enfa_calc_fun(locs = argo_loc,
species = "horridus",
reproductive_mode = "sexual",
mask_raster = ahor_bg_env)
#enfa for the asexual species
ahor_asexual_enfa <- enfa_calc_fun(locs = argo_loc,
species = "horridus",
reproductive_mode = "asexual",
mask_raster = ahor_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = ahor_sexual_enfa$m,
asex_marg = ahor_asexual_enfa$m,
full_species_name = "Argosarchus horridus")A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.
### 1) ENFA scatterplot
#access the relevant values for plotting
ahor_asexual_df <- ahor_asexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = ahor_asexual_enfa$pr)
ahor_sexual_df <- ahor_sexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = ahor_sexual_enfa$pr)
#asexual
enfa_hex_plot(ahor_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")
### 2) ENFA histogram
#asexual
hist(ahor_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
### 3) PCA of total e-space
total_climate_pca_plot(bg_env = ahor_bg_env, locs = argo_loc, genus = "Argosarchus", species = "horridus")binding character and factor vector, coercing into character vector
#pca
summary.list.aste <- species_pca_fun(loc.clim, "asteliaphasma")
#plot
aste_plot <- plot_clim_pca(summary.list.aste$loc.clim, summary.list.aste$summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. It's an annoying bug that hasn't been fixed yet.
#htmlwidgets::saveWidget(aste_plot, paste0(getwd(), "/plots/repro_mode_plots/asteliaphasma_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
aste_loc <- loc %>%
filter(genus == "asteliaphasma")
#use sourced plot_locs_leaflet script to plot localities
aste_map <- plot_locs_leaflet(aste_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(aste_map, url = paste0(getwd(), "/plots/repro_mode_plots/aste_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/aste_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7
Standard deviation 3.104 2.4714 1.27720 0.87566 0.76781 0.38573 0.19883
Proportion of Variance 0.507 0.3215 0.08586 0.04036 0.03103 0.00783 0.00208
Cumulative Proportion 0.507 0.8284 0.91427 0.95463 0.98566 0.99349 0.99557
PC8 PC9 PC10 PC11 PC12 PC13
Standard deviation 0.16377 0.1382 0.10580 0.09487 0.09187 0.06475
Proportion of Variance 0.00141 0.0010 0.00059 0.00047 0.00044 0.00022
Cumulative Proportion 0.99698 0.9980 0.99858 0.99905 0.99949 0.99971
PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.04986 0.03308 0.03192 0.02037 0.01726 0.01092
Proportion of Variance 0.00013 0.00006 0.00005 0.00002 0.00002 0.00001
Cumulative Proportion 0.99984 0.99990 0.99996 0.99998 0.99999 1.00000
loadings.aste <- summary.list.aste$summary.pca$rotation
knitr::kable(round(loadings.aste[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | -0.298 | -0.090 | -0.229 |
| chelsa_bioclims_NZ.2 | -0.288 | -0.043 | -0.320 |
| chelsa_bioclims_NZ.3 | -0.298 | -0.118 | -0.173 |
| chelsa_bioclims_NZ.4 | 0.198 | -0.312 | -0.110 |
| chelsa_bioclims_NZ.5 | 0.134 | -0.357 | -0.040 |
| chelsa_bioclims_NZ.6 | 0.209 | -0.281 | -0.137 |
| chelsa_bioclims_NZ.7 | -0.167 | -0.205 | 0.152 |
| chelsa_bioclims_NZ.8 | 0.143 | -0.352 | -0.071 |
| chelsa_bioclims_NZ.9 | 0.214 | -0.276 | -0.148 |
| chelsa_bioclims_NZ.10 | 0.187 | -0.284 | -0.161 |
| chelsa_bioclims_NZ.11 | 0.164 | -0.334 | -0.085 |
| chelsa_bioclims_NZ.12 | 0.250 | 0.211 | -0.238 |
| chelsa_bioclims_NZ.13 | 0.230 | 0.192 | -0.354 |
| chelsa_bioclims_NZ.14 | 0.266 | 0.202 | -0.053 |
| chelsa_bioclims_NZ.15 | -0.103 | 0.170 | -0.654 |
| chelsa_bioclims_NZ.16 | -0.295 | -0.154 | -0.044 |
| chelsa_bioclims_NZ.17 | 0.257 | 0.216 | -0.192 |
| chelsa_bioclims_NZ.18 | -0.297 | -0.125 | -0.157 |
| chelsa_bioclims_NZ.19 | -0.211 | -0.027 | -0.154 |
Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.
#get background env't for the species
ajuc_bg_env <- bg_env_crop(aste_loc,
species = "jucundum",
environment = w,
buffer = 0.5)no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf
#enfa for the sexual species
ajuc_sexual_enfa <- enfa_calc_fun(locs = aste_loc,
species = "jucundum",
reproductive_mode = "sexual",
mask_raster = ajuc_bg_env)
#enfa for the asexual species
ajuc_asexual_enfa <- enfa_calc_fun(locs = aste_loc,
species = "jucundum",
reproductive_mode = "asexual",
mask_raster = ajuc_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = ajuc_sexual_enfa$m,
asex_marg = ajuc_asexual_enfa$m,
full_species_name = "Asteliaphasma jucundum")A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.
### 1) ENFA scatterplot
#access the relevant values for plotting
ajuc_asexual_df <- ajuc_asexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = ajuc_asexual_enfa$pr)
ajuc_sexual_df <- ajuc_sexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = ajuc_sexual_enfa$pr)
#asexual
enfa_hex_plot(ajuc_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")
### 2) ENFA histogram
#asexual
hist(ajuc_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
### 3) PCA of total e-space
total_climate_pca_plot(bg_env = ajuc_bg_env, locs = aste_loc, genus = "Asteliophasma", species = "jucundum")binding character and factor vector, coercing into character vector
summary.list.clita <- species_pca_fun(loc.clim, "clitarchus")
clita_plot <- plot_clim_pca(summary.list.clita$loc.clim, summary.list.clita$summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. It's an annoying bug that hasn't been fixed yet.
#htmlwidgets::saveWidget(clita_plot, paste0(getwd(), "/plots/repro_mode_plots/clitarchus_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
clita_loc <- loc %>%
filter(genus == "clitarchus")
#use sourced plot_locs_leaflet script to plot localities
clita_map <- plot_locs_leaflet(clita_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(clita_map, url = paste0(getwd(), "/plots/repro_mode_plots/clita_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/clita_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7
Standard deviation 3.1671 2.2063 1.37300 1.10390 0.7268 0.47246 0.36559
Proportion of Variance 0.5279 0.2562 0.09922 0.06414 0.0278 0.01175 0.00703
Cumulative Proportion 0.5279 0.7841 0.88333 0.94747 0.9753 0.98702 0.99405
PC8 PC9 PC10 PC11 PC12 PC13
Standard deviation 0.22448 0.15341 0.12224 0.09218 0.07237 0.06626
Proportion of Variance 0.00265 0.00124 0.00079 0.00045 0.00028 0.00023
Cumulative Proportion 0.99670 0.99794 0.99873 0.99918 0.99945 0.99968
PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.04930 0.04095 0.03351 0.01839 0.01746 0.01160
Proportion of Variance 0.00013 0.00009 0.00006 0.00002 0.00002 0.00001
Cumulative Proportion 0.99981 0.99990 0.99996 0.99998 0.99999 1.00000
loadings.clita <- summary.list.clita$summary.pca$rotation
knitr::kable(round(loadings.clita[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | -0.277 | -0.151 | 0.244 |
| chelsa_bioclims_NZ.2 | -0.267 | -0.106 | 0.338 |
| chelsa_bioclims_NZ.3 | -0.281 | -0.172 | 0.173 |
| chelsa_bioclims_NZ.4 | -0.224 | 0.311 | -0.091 |
| chelsa_bioclims_NZ.5 | -0.245 | 0.256 | -0.157 |
| chelsa_bioclims_NZ.6 | -0.215 | 0.312 | -0.002 |
| chelsa_bioclims_NZ.7 | -0.133 | -0.109 | -0.294 |
| chelsa_bioclims_NZ.8 | -0.245 | 0.260 | -0.147 |
| chelsa_bioclims_NZ.9 | -0.213 | 0.316 | -0.013 |
| chelsa_bioclims_NZ.10 | -0.216 | 0.294 | -0.001 |
| chelsa_bioclims_NZ.11 | -0.233 | 0.274 | -0.139 |
| chelsa_bioclims_NZ.12 | 0.234 | 0.245 | 0.241 |
| chelsa_bioclims_NZ.13 | 0.231 | 0.232 | 0.203 |
| chelsa_bioclims_NZ.14 | 0.205 | 0.262 | 0.271 |
| chelsa_bioclims_NZ.15 | -0.182 | 0.018 | 0.577 |
| chelsa_bioclims_NZ.16 | -0.280 | -0.203 | 0.071 |
| chelsa_bioclims_NZ.17 | 0.228 | 0.259 | 0.266 |
| chelsa_bioclims_NZ.18 | -0.192 | -0.204 | 0.060 |
| chelsa_bioclims_NZ.19 | -0.204 | -0.043 | 0.227 |
Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.
#get background env't for the species
choo_bg_env <- bg_env_crop(clita_loc,
species = "hookeri",
environment = w,
buffer = 0.5)
#enfa for the sexual species
choo_sexual_enfa <- enfa_calc_fun(locs = clita_loc,
species = "hookeri",
reproductive_mode = "sexual",
mask_raster = choo_bg_env)
#enfa for the asexual species
choo_asexual_enfa <- enfa_calc_fun(locs = clita_loc,
species = "hookeri",
reproductive_mode = "asexual",
mask_raster = choo_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = choo_sexual_enfa$m,
asex_marg = choo_asexual_enfa$m,
full_species_name = "Clitarchus hookeri")A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.
### 1) ENFA scatterplot
#access the relevant values for plotting
choo_asexual_df <- choo_asexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = choo_asexual_enfa$pr)
choo_sexual_df <- choo_sexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = choo_sexual_enfa$pr)
#asexual
enfa_hex_plot(choo_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")
### 2) ENFA histogram
#asexual
hist(choo_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
### 3) PCA of total e-space
total_climate_pca_plot(bg_env = choo_bg_env, locs = aste_loc, genus = "Clitarchus", species = "hookeri")binding character and factor vector, coercing into character vector
summary.list.micra <- species_pca_fun(loc.clim, "micrarchus")
micra_plot <- plot_clim_pca(summary.list.micra$loc.clim, summary.list.micra$summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. It's an annoying bug that hasn't been fixed yet.
#htmlwidgets::saveWidget(micra_plot, paste0(getwd(), "/plots/repro_mode_plots/micrarchus_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
micra_loc <- loc %>%
filter(genus == "micrarchus")
#use sourced plot_locs_leaflet script to plot localities
micra_map <- plot_locs_leaflet(micra_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(micra_map, url = paste0(getwd(), "/plots/repro_mode_plots/micra_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/micra_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7
Standard deviation 3.4431 2.0743 1.31392 0.94676 0.29817 0.2263 0.20567
Proportion of Variance 0.6239 0.2265 0.09086 0.04718 0.00468 0.0027 0.00223
Cumulative Proportion 0.6239 0.8504 0.94126 0.98843 0.99311 0.9958 0.99803
PC8 PC9 PC10 PC11 PC12 PC13
Standard deviation 0.13292 0.09179 0.08909 0.03357 0.02983 0.02512
Proportion of Variance 0.00093 0.00044 0.00042 0.00006 0.00005 0.00003
Cumulative Proportion 0.99896 0.99941 0.99983 0.99988 0.99993 0.99996
PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.01449 0.01179 0.01068 0.009585 0.008206 0.007186
Proportion of Variance 0.00001 0.00001 0.00001 0.000000 0.000000 0.000000
Cumulative Proportion 0.99998 0.99998 0.99999 0.999990 1.000000 1.000000
loadings.micra <- summary.list.micra$summary.pca$rotation
knitr::kable(round(loadings.micra[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | 0.268 | 0.116 | -0.222 |
| chelsa_bioclims_NZ.2 | 0.266 | 0.086 | -0.270 |
| chelsa_bioclims_NZ.3 | 0.268 | 0.143 | -0.178 |
| chelsa_bioclims_NZ.4 | -0.242 | 0.257 | -0.101 |
| chelsa_bioclims_NZ.5 | -0.237 | 0.254 | -0.107 |
| chelsa_bioclims_NZ.6 | -0.239 | 0.251 | -0.153 |
| chelsa_bioclims_NZ.7 | 0.145 | -0.116 | -0.015 |
| chelsa_bioclims_NZ.8 | -0.237 | 0.254 | -0.118 |
| chelsa_bioclims_NZ.9 | -0.239 | 0.253 | -0.144 |
| chelsa_bioclims_NZ.10 | -0.240 | 0.250 | -0.147 |
| chelsa_bioclims_NZ.11 | -0.229 | 0.252 | -0.102 |
| chelsa_bioclims_NZ.12 | -0.167 | -0.323 | -0.354 |
| chelsa_bioclims_NZ.13 | -0.140 | -0.310 | -0.429 |
| chelsa_bioclims_NZ.14 | -0.195 | -0.325 | -0.202 |
| chelsa_bioclims_NZ.15 | 0.247 | 0.004 | -0.397 |
| chelsa_bioclims_NZ.16 | 0.268 | 0.170 | -0.110 |
| chelsa_bioclims_NZ.17 | -0.190 | -0.311 | -0.294 |
| chelsa_bioclims_NZ.18 | 0.210 | 0.226 | -0.306 |
| chelsa_bioclims_NZ.19 | 0.267 | 0.119 | -0.185 |
summary.list.nive <- species_pca_fun(loc.clim, "niveaphasma")
nive_plot <- plot_clim_pca(summary.list.nive$loc.clim, summary.list.nive$summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. It's an annoying bug that hasn't been fixed yet.
#htmlwidgets::saveWidget(nive_plot, paste0(getwd(), "/plots/repro_mode_plots/niveaphasma_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
nive_loc <- loc %>%
filter(genus == "niveaphasma")
#use sourced plot_locs_leaflet script to plot localities
nive_map <- plot_locs_leaflet(nive_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(nive_map, url = paste0(getwd(), "/plots/repro_mode_plots/nive_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/nive_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7
Standard deviation 2.9922 2.4493 1.5457 0.94013 0.73205 0.39659 0.24024
Proportion of Variance 0.4712 0.3157 0.1258 0.04652 0.02821 0.00828 0.00304
Cumulative Proportion 0.4712 0.7870 0.9127 0.95923 0.98743 0.99571 0.99875
PC8 PC9 PC10 PC11 PC12 PC13
Standard deviation 0.11345 0.06508 0.04982 0.04382 0.02749 0.02222
Proportion of Variance 0.00068 0.00022 0.00013 0.00010 0.00004 0.00003
Cumulative Proportion 0.99943 0.99965 0.99978 0.99988 0.99992 0.99995
PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.01856 0.01631 0.01216 0.01109 0.008331 0.006935
Proportion of Variance 0.00002 0.00001 0.00001 0.00001 0.000000 0.000000
Cumulative Proportion 0.99997 0.99998 0.99999 0.99999 1.000000 1.000000
loadings.nive <- summary.list.nive$summary.pca$rotation
knitr::kable(round(loadings.nive[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | 0.091 | 0.374 | 0.184 |
| chelsa_bioclims_NZ.2 | -0.002 | 0.341 | 0.347 |
| chelsa_bioclims_NZ.3 | 0.156 | 0.358 | 0.053 |
| chelsa_bioclims_NZ.4 | 0.297 | -0.154 | 0.159 |
| chelsa_bioclims_NZ.5 | 0.294 | -0.151 | 0.172 |
| chelsa_bioclims_NZ.6 | 0.297 | -0.158 | 0.135 |
| chelsa_bioclims_NZ.7 | -0.196 | 0.021 | 0.110 |
| chelsa_bioclims_NZ.8 | 0.296 | -0.143 | 0.177 |
| chelsa_bioclims_NZ.9 | 0.298 | -0.158 | 0.146 |
| chelsa_bioclims_NZ.10 | 0.300 | -0.143 | 0.132 |
| chelsa_bioclims_NZ.11 | 0.313 | -0.119 | 0.116 |
| chelsa_bioclims_NZ.12 | -0.234 | -0.233 | 0.274 |
| chelsa_bioclims_NZ.13 | -0.216 | -0.203 | 0.338 |
| chelsa_bioclims_NZ.14 | -0.241 | -0.238 | 0.233 |
| chelsa_bioclims_NZ.15 | -0.119 | 0.211 | 0.496 |
| chelsa_bioclims_NZ.16 | 0.176 | 0.346 | -0.012 |
| chelsa_bioclims_NZ.17 | -0.234 | -0.243 | 0.252 |
| chelsa_bioclims_NZ.18 | -0.171 | 0.256 | 0.192 |
| chelsa_bioclims_NZ.19 | 0.116 | 0.151 | 0.289 |
Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.
#get background env't for the species
nive_bg_env <- bg_env_crop(nive_loc,
species = "annulata",
environment = w,
buffer = 0.5)
#enfa for the sexual species
nive_sexual_enfa <- enfa_calc_fun(locs = nive_loc,
species = "annulata",
reproductive_mode = "sexual",
mask_raster = nive_bg_env)
#enfa for the asexual species
nive_asexual_enfa <- enfa_calc_fun(locs = nive_loc,
species = "annulata",
reproductive_mode = "asexual",
mask_raster = nive_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = nive_sexual_enfa$m,
asex_marg = nive_asexual_enfa$m,
full_species_name = "Niveaphasma annulata")A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.
### 1) ENFA scatterplot
#access the relevant values for plotting
nive_asexual_df <- nive_asexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = nive_asexual_enfa$pr)
nive_sexual_df <- nive_sexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = nive_sexual_enfa$pr)
#asexual
enfa_hex_plot(nive_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")
### 2) ENFA histogram
#asexual
hist(nive_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
### 3) PCA of total e-space
total_climate_pca_plot(bg_env = nive_bg_env, locs = nive_loc, genus = "Niveaphasma", species = "annulata")binding character and factor vector, coercing into character vector
summary.list.spin <- species_pca_fun(loc.clim, "spinotectarchus")
spin_plot <- plot_clim_pca(summary.list.spin$loc.clim, summary.list.spin$summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. It's an annoying bug that hasn't been fixed yet.
#htmlwidgets::saveWidget(spin_plot, paste0(getwd(), "/plots/repro_mode_plots/spinotectarchus_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
spin_loc <- loc %>%
filter(genus == "spinotectarchus")
#use sourced plot_locs_leaflet script to plot localities
spin_map <- plot_locs_leaflet(spin_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(spin_map, url = paste0(getwd(), "/plots/repro_mode_plots/spin_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/spin_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7
Standard deviation 3.2696 2.2078 1.34284 0.86560 0.74380 0.4655 0.17853
Proportion of Variance 0.5626 0.2566 0.09491 0.03944 0.02912 0.0114 0.00168
Cumulative Proportion 0.5626 0.8192 0.91410 0.95354 0.98265 0.9941 0.99573
PC8 PC9 PC10 PC11 PC12 PC13
Standard deviation 0.16747 0.12847 0.11756 0.10182 0.07402 0.05291
Proportion of Variance 0.00148 0.00087 0.00073 0.00055 0.00029 0.00015
Cumulative Proportion 0.99721 0.99808 0.99881 0.99935 0.99964 0.99979
PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.04093 0.03100 0.02372 0.02118 0.01686 0.01123
Proportion of Variance 0.00009 0.00005 0.00003 0.00002 0.00001 0.00001
Cumulative Proportion 0.99987 0.99993 0.99995 0.99998 0.99999 1.00000
loadings.spin <- summary.list.spin$summary.pca$rotation
knitr::kable(round(loadings.spin[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | -0.279 | -0.110 | -0.223 |
| chelsa_bioclims_NZ.2 | -0.273 | -0.057 | -0.289 |
| chelsa_bioclims_NZ.3 | -0.277 | -0.148 | -0.178 |
| chelsa_bioclims_NZ.4 | 0.243 | -0.262 | -0.126 |
| chelsa_bioclims_NZ.5 | 0.199 | -0.325 | -0.112 |
| chelsa_bioclims_NZ.6 | 0.246 | -0.230 | -0.070 |
| chelsa_bioclims_NZ.7 | -0.133 | -0.273 | -0.032 |
| chelsa_bioclims_NZ.8 | 0.209 | -0.314 | -0.127 |
| chelsa_bioclims_NZ.9 | 0.254 | -0.217 | -0.106 |
| chelsa_bioclims_NZ.10 | 0.231 | -0.236 | -0.144 |
| chelsa_bioclims_NZ.11 | 0.215 | -0.303 | -0.135 |
| chelsa_bioclims_NZ.12 | 0.199 | 0.273 | -0.319 |
| chelsa_bioclims_NZ.13 | 0.174 | 0.242 | -0.434 |
| chelsa_bioclims_NZ.14 | 0.225 | 0.267 | -0.030 |
| chelsa_bioclims_NZ.15 | -0.165 | 0.138 | -0.579 |
| chelsa_bioclims_NZ.16 | -0.274 | -0.193 | -0.055 |
| chelsa_bioclims_NZ.17 | 0.217 | 0.281 | -0.234 |
| chelsa_bioclims_NZ.18 | -0.272 | -0.164 | -0.178 |
| chelsa_bioclims_NZ.19 | -0.206 | -0.046 | -0.140 |
summary.list.tect <- species_pca_fun(loc.clim, "tectarchus")
tect_plot <- plot_clim_pca(summary.list.tect$loc.clim, summary.list.tect$summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. It's an annoying bug that hasn't been fixed yet.
#htmlwidgets::saveWidget(tect_plot, paste0(getwd(), "/plots/repro_mode_plots/tectarchus_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
tect_loc <- loc %>%
filter(genus == "tectarchus")
#use sourced plot_locs_leaflet script to plot localities
tect_map <- plot_locs_leaflet(tect_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(tect_map, url = paste0(getwd(), "/plots/repro_mode_plots/tect_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/tect_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7
Standard deviation 2.8210 2.5950 1.6213 0.97360 0.59328 0.41536 0.34340
Proportion of Variance 0.4188 0.3544 0.1383 0.04989 0.01853 0.00908 0.00621
Cumulative Proportion 0.4188 0.7733 0.9116 0.96150 0.98002 0.98910 0.99531
PC8 PC9 PC10 PC11 PC12 PC13
Standard deviation 0.23982 0.11480 0.08517 0.06494 0.04977 0.03772
Proportion of Variance 0.00303 0.00069 0.00038 0.00022 0.00013 0.00007
Cumulative Proportion 0.99833 0.99903 0.99941 0.99963 0.99976 0.99984
PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.03413 0.02872 0.02484 0.01571 0.01133 0.01044
Proportion of Variance 0.00006 0.00004 0.00003 0.00001 0.00001 0.00001
Cumulative Proportion 0.99990 0.99994 0.99997 0.99999 0.99999 1.00000
loadings.tect <- summary.list.tect$summary.pca$rotation
knitr::kable(round(loadings.tect[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | 0.333 | 0.029 | -0.197 |
| chelsa_bioclims_NZ.2 | 0.310 | 0.019 | -0.288 |
| chelsa_bioclims_NZ.3 | 0.345 | 0.030 | -0.114 |
| chelsa_bioclims_NZ.4 | 0.048 | -0.380 | 0.021 |
| chelsa_bioclims_NZ.5 | 0.070 | -0.367 | 0.089 |
| chelsa_bioclims_NZ.6 | 0.030 | -0.379 | -0.011 |
| chelsa_bioclims_NZ.7 | 0.102 | 0.106 | 0.275 |
| chelsa_bioclims_NZ.8 | 0.071 | -0.369 | 0.078 |
| chelsa_bioclims_NZ.9 | 0.027 | -0.379 | -0.010 |
| chelsa_bioclims_NZ.10 | 0.018 | -0.379 | -0.015 |
| chelsa_bioclims_NZ.11 | 0.082 | -0.352 | 0.081 |
| chelsa_bioclims_NZ.12 | -0.279 | -0.060 | -0.352 |
| chelsa_bioclims_NZ.13 | -0.252 | -0.058 | -0.362 |
| chelsa_bioclims_NZ.14 | -0.300 | -0.045 | -0.285 |
| chelsa_bioclims_NZ.15 | 0.215 | -0.011 | -0.481 |
| chelsa_bioclims_NZ.16 | 0.351 | 0.042 | -0.021 |
| chelsa_bioclims_NZ.17 | -0.285 | -0.062 | -0.343 |
| chelsa_bioclims_NZ.18 | 0.273 | 0.008 | -0.209 |
| chelsa_bioclims_NZ.19 | 0.293 | 0.046 | -0.197 |
Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.
This is for Tectarchus ovobessus.
#get background env't for the species
tect_ovo_bg_env <- bg_env_crop(tect_loc,
species = "ovobessus",
environment = w,
buffer = 0.5)
#enfa for the sexual species
tect_ovo_sexual_enfa <- enfa_calc_fun(locs = tect_loc,
species = "ovobessus",
reproductive_mode = "sexual",
mask_raster = tect_ovo_bg_env)
#enfa for the asexual species
tect_ovo_asexual_enfa <- enfa_calc_fun(locs = tect_loc,
species = "ovobessus",
reproductive_mode = "asexual",
mask_raster = tect_ovo_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = tect_ovo_sexual_enfa$m,
asex_marg = tect_ovo_asexual_enfa$m,
full_species_name = "Tectarchus ovobessus")A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.
### 1) ENFA scatterplot
#access the relevant values for plotting
tect_ovo_asexual_df <- tect_ovo_asexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = tect_ovo_asexual_enfa$pr)
tect_ovo_sexual_df <- tect_ovo_sexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = tect_ovo_sexual_enfa$pr)
#asexual
enfa_hex_plot(tect_ovo_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")
### 2) ENFA histogram
#asexual
hist(tect_ovo_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
### 3) PCA of total e-space
total_climate_pca_plot(bg_env = tect_ovo_bg_env, locs = tect_loc, genus = "Tectarchus", species = "ovobessus")binding character and factor vector, coercing into character vector
This is an enfa for Tectarchus huttoni.
###Only need to get bg env't if you're not running the previous chunk
#get background env't for the species
tect_hutt_bg_env <- bg_env_crop(tect_loc,
species = "huttoni",
environment = w,
buffer = 0.5)
#enfa for the sexual species
tect_hutt_sexual_enfa <- enfa_calc_fun(locs = tect_loc,
species = "huttoni",
reproductive_mode = "sexual",
mask_raster = tect_hutt_bg_env)
#enfa for the asexual species
tect_hutt_asexual_enfa <- enfa_calc_fun(locs = tect_loc,
species = "huttoni",
reproductive_mode = "asexual",
mask_raster = tect_hutt_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = tect_hutt_sexual_enfa$m,
asex_marg = tect_hutt_asexual_enfa$m,
full_species_name = "Tectarchus huttoni")A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.
### 1) ENFA scatterplot
#access the relevant values for plotting
tect_hutt_asexual_df <- tect_hutt_asexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = tect_hutt_asexual_enfa$pr)
tect_hutt_sexual_df <- tect_hutt_sexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = tect_hutt_sexual_enfa$pr)
#asexual
enfa_hex_plot(tect_hutt_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")
### 2) ENFA histogram
#asexual
hist(tect_hutt_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
### 3) PCA of total e-space
total_climate_pca_plot(bg_env = tect_hutt_bg_env, locs = tect_loc, genus = "Tectarchus", species = "huttoni")binding character and factor vector, coercing into character vector
Nothing. Only one locality.